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Abstract 

The u» of |n&it bind o ptimis a tion iljnitkai is aww don pi is wOl established as 
a practical approach to aerodynamic design. A typical procedure uaaa a emulation scheme to 
evaluate the objective f un cti on (from tbs eppr estimate states) and ita padieat, thus passes this 
information to an optimisation algorithm. Once the rimulation scheme (CFD few solver) has 
been selected ami used to provide approtimata function evaluations, there are tevaral possible 
approaches to the problem of computing gradients. One popular method is to dffiaentiate 
the simulation scheme and compute design sensitivities that an then used to obtain gradients. 
Although this black-box approach has many advantages in shape optimisation problems, me 
must compute mesh eentitivhies in order to compute the derip senritivity. 

In this paper, we pnweat an alternative approach using the PDE sensitivity equation to 
develop algorithms for computing gradients. This approach hss the advantage that met h eenri- 
tivitiee need not be computed. Moreover, whan it is possible to use the CPD scheme for both 
the forward problem and the sensitivity equation, then there am computational advantages. An 
apparent disadvantage of this approach is that it does not always produce consistent deriva- 
tives. However, for a proper c ombin ation of discretisation schemes, me can show asymptotic 
consistency seder ma rt xednemcat, which is often sufRcient to psraatee convergence of the 
optimal design algorithm. In particular, we show that when asymptotically nontisteut schemes 
an combined with a trust-region optimisation algorithm, the resulting optimal design method 
converges. We denote this approach as the Miwitisity eguetion method. 

The sensitivity equation method is presented, convergence results are given and the approach 
h illustrated on two optimal design problems involving shocks. 
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1 Introduction 


Optimal design problems consist of selecting design parameters for a system in order 
to optimize a given design objective, usually constrained to satisfy a partial differential 
equation. In many of these problems, design parameters describe the shape of an 
object. Examples of these shape optimization problems include drag reduction [21], 
[22], weight minimization [14], optimal sensor /actuator placement [6], airfoil design 
[16], [17], [18], [19] and the design of wind tunnel elements [15]. 

Traditionally, approximate solutions of these problems are found by “cut and try” 
methods, combining a designer’s engineering experience with repeated experimental 
testing. This is often expensive, motivating computational methods which compute 
the optimal design directly. These methods require defining an objective function 
and an appropriate PDE model of the states of the system. A comparison of several 
optimal design methods may be found in [13]. 

Many popular approaches couple a gradient-based optimization algorithm with 
function evaluations provided by a proven simulation scheme. One of the disadvan- 
tages of these approaches is the expense of computing the gradient. Using finite 
differences is often too costly, even if appropriate step sizes can be found and the sim- 
ulation scheme can take advantage of “nearby’’ solutions (as is the case with iterative 
solvers for nonlinear equations). 

Two strategies for alleviating the computational expense of gradient evaluations 
are adjoint variables [17] and design sensitivities [14]. Adjoint methods are advanta- 
geous when either the problem is self-adjoint or there are a large number of design 
parameters. However, when there are relatively few design parameters, wring design 
sensitivities, quantities which describe the influence of the design parameters on the 
states of the system, is an attractive alternative. In addition to efficient gradient 
computations, they can be used in some problems to construct an effective update of 
the approximate Hessian for quasi-Newton optimization algorithms, e.g. [10]. 

A standard approach often used to compute the design sensitivities is based on 
(implicitly) differentiating the simulation scheme (for the states) with respect to the 
design variables. Using the chain rule to carry out this calculation, results in an 
efficient numerical scheme for the sensitivities. The efficiency arises from reusing 
many of the quantities computed in the simulation scheme. In fact, the “inversion” 
of the system matrix (i.e. the matrix factorization) can often be reused. 
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A disadvantage of this approach is that for shape optimization prob lems , the dis- 
cretization is parameter dependent. Thus, derivatives of the discretization (mesh 
sensitivities) are required for each shape parameter. Depending on the simulation 
scheme used for the states, determining the discretization can require the solution of 
a partial differentia] equation (as is the case for finite difference solutions of viscous 
flow problems (26]). This requires a strategy for computing the mesh sensitivities 
(20], or for computing an approximation to them (24], (25]. 

Another approach to finding design sensitivities relies on approximating the par- 
tial differential equation, known as the sensitivity equation. This equation is obtained 
by implicitly differentiating the (infinite dimenskn.al) state equation with respect to 
each design parameter. As shown in (2], using the same numerical scheme to ap- 
proximate the sensitivity equation which is used to approximate the states, leads to 
an efficient scheme with similar computational advantages as the design sensitivity 
approach described above. Furthermore, since the discretization is applied directly 
to the sensitivity equations, no sensitivity of the mesh is required. The sensitivity 
equation is always linear in the design sensitivity, even if the state equation is non- 
linear. Since there is no requirement to use the same numerical scheme, it is possible 
to gain additional computational savings by using a scheme which tak es advantage 
of the linearity in the sensitivity equations. 

An apparent disadvantage of this approach is that it does not compute consistent 
derivatives. In other words, the sensitivity equation approach does not capture the 
sensitivity of the truncation errors in the scheme. Thus, there is a concern that 
providing an optimization algorithm with an approximation of the gradient of the 
infinite dimensional objective function instead of the gradient of the approximate 
objective function would cause the algorithm to fail. One might expect, however, 
that if the gradients are “close enough” to the true gradients, then the optimization 
algorithm should still converge. We show that this convergence can be established if 
one combines compatible simulation and optimization schemes. 

TV ust- region optimization algorithms are constructed to be globally convergent by 
minimizing a model of the objective function in a region where the model is “trusted” . 
This leads to robust algorithms capable of handling inaccuracies in the model. In fact, 
convergence results have been given for these algorithms when the model is based on 
inaccurate gradient information (7], (8). The results hold provided the grad'ents satisfy 
a given error condition. Therefore, it is natural to consider an optimal design method 
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which couples a trust-region optimization algorithm with gradients computed using 
the sensitivity equation approach. We denote this combined sensitivity /trust-region 
algorithm by the sensitivity equation method (SEM). 

In this work, we present and analyze the sensitivity equation method. The method 
can be applied to a wide class of optimal design problems, including those mentioned 
above, however, we focus mi the particular example of shape optimization of Euler 
flows in order to illustrate the method. In the next section, we describe two design 
problems. In Section 3, we present the sensitivity equation method including the 
trust-region algorithm and the use of the sensitivity equation to find the design sensi- 
tivities. Furthermore, we compare various numerical approximations of the sensitivity 
equation w.*h approaches based on the discretized equations. Section 4 discusses a 
number of convergence issues and includes a convergence theorem for the sensitivity 
equation method. In Section 5, we use a one dimensional duct design problem to de- 
scribe the implementation of the sensitivity equation approach. Finally, we describe 
the implementation and perform shape optimization for a two dimensional forebody 
simulator design problem where the steady-state Euler equations are used to model 
the state variables. 

2 Illustrative Examples 

We present two optimal design problems below which are used to illustrate the sen- 
sitivity equation method. These problems consist of determining shape parameters 
which produce a solution to the Euler equations that matches a desired flow “as 
closely as possible.” The first problem is motivated by the design of a wind tunnel 
element in order to produce a desired flow in the test section. We study a two di- 
mensional analogue of this problem. The second problem consists of prescribing the 
cross-sectional area of a one dimensional duct to produce a duct flow which matches 
a desired flow profile. This problem was used by Frank and Shubin (13) in their study 
of optimal design. 

2.1 Forebody Simulator Design Problem 

The Arnold Engineering Development Center (AEDC) operates a free-jet test facility 
which is used for full-scale testing of commercial and military aircraft engines. Engines 
arc evaluated for performance and safety under various free flight conditions. While 
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Figure 1: Forebody Simulator Design Problem 

this facility is large enough to house engines, it is not large enough to house an entire 
aircraft forebody. Thus, the effect of the aircraft forebody on the engine inlet flow 
profile must be simulated. One way of doing this is to replace the actual forebody by 
a smaller object, called a forebody simulator (FBS). The use of the FBS is illustrated 
in Figure 1. The FBS design problem is to specify the shape of this FBS so that 
it produces an engine inlet flow profile which is as close to some desired profile as 
possible [15 j. The desired profile can be determined by performing either a wind 
tunnel simulation or a computational simulation of a model configuration resembling 
a test condition of the aircraft engine. 

In order to demonstrate the applicability of the SEM, we consider a two dimen* 
sional analogue of this problem. This problem, depicted in Figure 2, is to find the 
shape of the curve I\ which produces an outflow that matches the outflow generated 
by the original (longer) forebody as closely as possible. The flow, Q (consisting of the 
density p, the momentum put -f pvj and the sum of the internal and kinetic energy 
E) is modeled using the steady state Euler equations, 
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Shortened Forebody (Shape to be Determined) 


Figure 2: 2D Fore body Simulator Design Problem 
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The pressure P it related to the elements of Q by 


( 1 ) 

( 2 ) 


P*(7-l)[s-^p(u 3 + t; 3 )], (3) 

where 7 is the ratio of specific heats (7 ■= 1.4 for air). Given a forebody simulator 
shape T. the flow <?(r) is determined by solving the Euler equations ( 1 ) in the test 
cell domain ft(T) subject to the boundary conditions (for supersonic flow): 


Q 

(u, v) • n 



Qim at the test cell inflow, 

0 and 

0 at the walls (no flow penetration), 
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(4) 

(5) 

( 6 ) 



where n and t axe the Doimal and tangential vectors at the boundary, respectively. 
The set of admissible forebody simulator shapes is 

a - {r e c‘M)| r(a) « r., r (*) = r b and r(x) > r„ Vx € (a,*)} . (7) 

A statement of the design problem is given below. 

Problem 2.1 (Forebody Simulator Design) Let C) be a desired flow at the out- 
flow (engine inlet), 

S= {{*,y)\* = b, T\,<y<c}. (8) 

Define the objective function 

j(r>- j)Q(T)-<S\l ds, ( 9 ) 

where <?(T) represents the solution of (1) with boundary conditions (4)- (6) in the test 
cell fl(T). The forebody simulator design problem is to find T. € A such that 

J(T.)<J(T) for all T € -4. (10) 

Closed form solutions to (1) with (4)-(6) are available only for special domains. 
Therefore, we consider approximate solutions of (1) and hence the approximation of 
Problem 2.1. 

The discretization is performed by selecting mesh points in the flow domain fi(P) 
where the flow variables will be approximated. It is desirable to select this mesh 
in such a way that the points are more dense in regions where flow gradients are 
expected to be “large” (in order to have more accurate differencing) and more coarse 
in regions where the flow is nearly constant (in order to save computer time). Other 
issues such as selecting points with no sharp changes in density and with sufficient 
resolution to treat the boundary conditions, make the mesh generation a science in 
and of itself (see e.g. [26]). 

Another constraint on the discretization, to simplify the implementation of a finite 
difference scheme, is to use a regular mesh, i.e. a mesh where there exists a bijective 
map taking the mesh points to a lattice of points in the computational space. For 
example, suppose that M is a C l mapping, 

M (x,y) -> (Z,r\), (11) 

then derivatives in the physical space are easily approximated on the lattice using 
the chain rule. Denoting the Jacobian of the mapping by Jjh, the transformed Euler 
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(13) 


Q = JmQ, = „d v«g„ + |„. (u) 

A standard finite difference scheme, developed by Beam and Warming [1] is used 
to approximate the transformed equations. The scheme introduces a time variable, t 
as a means of iterating an initial guess for the solution, to a solution of the steady 
state equations. Second and fourth order artificial dissipation terms are added for 
stability, represented by and respectively. This scheme is implemented 
in the PARC2D code {9]. Several implementation issues are discussed briefly below 
which are referred to in later sections. Readers interested in more code details or the 
actual expressions used for 4^*1 and should consult [9]. 

The difference scheme produces a system of equations for the update of the flow 
variables, AQ n . Thus, the solution at the nth iteration, Q n is determined from 


<J" * Q'- 1 + A#"” 1 . 


(15) 


The system matrix produced by the approximation above is quite large due to differ- 
encing in each direction. However, this problem is circumvented us ing an approximate 
factorization into a product of two matrices, each corresponding to differencing in one 
of the lattice directions. The final system has the form: 

'/ + A t6 ( A* - V< > + *[ 4) ) A iJ M \ x 
/ + AtSiB" - V, (tfW + *<«>) A^Jm] A$ n * 

- A t6 ( P* - A t6,Gr 

+ A tV ( (¥<’> - tJ%V c ) A ( (J m Q») 

+ AtV v (*»-¥*A v V,)A v (j M Q«) 1 (16) 

where 

“ d ( 1? ) 
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The subscripted terms 6 , V and A represent the central, backward and forward 
difference operators, respectively, in the lattice direction indicated by the subscript. 
The converged solution is denoted by Q*(x,y) = Q n * (Af(x, y)). 

We introduce Bezier curves to parameterize the forebody simulator. Bezier poly- 
nomials possess several nice properties when used in approximations. The most im- 
portant for the examples presented here are the ecnvtz hull and endpoint interpolation 
properties (see e.g. Farm [12]). For this problem, we consider a set of two parameter, 
q = (?\g 2 ), Bezier curves 

5 = {rec'[o,i]| n*)- (r.W.r, (.;«)), r,(, !t )sr., .«(o,i], ? €i>} (is) 

where 

T,(s) * oBb^(s) + 0.6B M (s)+0.8B 3t3 (s) + 6B ai3 (s), (19) 

r„(s;g) = r.B 0l3 (s) + q'Bufs) + g 2 Bj,j(a) + rhB 3>3 (s), ( 20 ) 

and 

£«>(*) = ^ j *’(1 - x ) r -\ ( 21 ) 

We also assume a as 0.5 and b * 1.0. We can now introduce the approximate forebody 
simulator design problem. 

Problem 2.2 (Approximate Forebody Simulator Design) Let {&}' be de- 
sired flow measurements at S. We assume that the data measurements' are given 
at the quadrature points, otherwise interpolation must be used. Define the objective 

function 

J^ r ) = X>|Q''(*«r)-<9 ( £, ( 22 ) 

^•1 

where Q N (zi;T) represents the approximate solution to (1) in the domain fl(r) at the 
quadrature point x,. The approximate forebody simulator design problem is to find 
T. € B such that 

4/V.) < J?(T) for all reS. (23) 

Let 

a -{(?’,«’)€ R’ | r(-,,V)€B}, (24) 

then the problem can be equivalently stated as finding (ql,q*) € Q such that 

J?(<llll)<J, N {q\q>) for all (o\g 7 ) € Q. (25) 
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2.2 Duct Design Problem 


This problem consists of designing the cross-sectional area of a one-dimensional duct 
such that, under specified inlet and outlet conditions, it produces a flow which is as 
close to a desired transonic flow as possible. The governing conservation laws (steady 
state continuity, momentum and energy equations) can be reduced to a single two- 
point boundary value problem (BVP) for the velocity. It was shown in (13) that the 
velocity u, is the solution of 

J“/0*)+$(«.A) = 0, (26) 

u( 0 ) = uja and u(l) *= tw, 

where uin and tteut are the velocities at the inlet and outlet of the duct, A is the 
cross-sectional area of the duct, 

/(„) = » + £, ,(u, ,4) = I M< 1 f.lzl, (27) 

where fl and 7 axe flow constants taken to he 1.14 and 1 . 4 , respectively. The Rankine- 
Hngoniot condition yields the speed of sound as u, = y/H. Unique solutions of this 
BVP are guaranteed for monotone area functions, therefore, cross-sectional areas, A, 
are restricted to 

A - {a € C^O, 1) A(0) = A*, A(l) a A wi and ~A(x) > 0, Vx € (0, 1) J (28) 

for fixed inlet and outlet areas of Au, and A^t* We now describe the optimal design 
problem. 

Problem 2.3 (Duct Design) Let u(-) € L 3 ( 0 , 1 ) be a desired transonic flow profile 
for the duct and define the objective function by 

J[A)m l (u(x;A)-u(x )] 8 dx (29) 

where ti ( •; A) is the solution to (26) corresponding to A. The optimal design problem 
is to find on A, € A such that 

J{A.)<J(A) for all A € A. (30) 

While the BVP has a closed form solution (13), we consider approximations of 
(26) and consequently of Problem 2.3 in order to study the more general case. We 
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begin by discretizing the duct length into N cells (of length h = jj) with centers, 
x j = 0 “ j ** 1 1 • • • ♦ AT nod define to be the average velocity in the jth cell, 
i.e. 

u i (•*) * | J f m **i *(*> A )dx. (31) 

f + jst 

A system of nonlinear equations for u w (A) = {uf (A)| i can be found by integrating 
(26) over each cell. 


l (“ (f> ± |m)) - / (" (*j - M)) 


+ J «(A), A(x,)) = 0, ; = 


.,N, 
(32) 

where it was assumed that was nearly constant over each cell. An approxima- 
tion to u N is found by replacing the fluxes at the cell edges, / (u (xj + *)) , using the 
cell center values /, * /(uf ) and £+i = /(u£, ). Two standard first order “Godunov 
type” methods are the Enquist-Osher scheme 

fi - fi 


(“M)) 


** Ff+l/2 


fi 

/(«.) 

fi + fj+i - /(«,) 




,n w jy 

l s ’ u i+i 


> 


< u 4 < u&,; 

u ?+i <«*<«?• 


(33) 


and the artificial viscosity scheme 

f (u (xj -h * F$ /a = - 2 (/ j+1 + fi -q(u j+ , - Uj )) , 


(34) 


where a has been selected as 1 for this study. These approximations were used in 
[13], but are included above for completeness. 

We turn now to the approximation of the cross-sectional area A. The space A 
is replaced by a subset of Bezier quadratic polynomials. The properties of Bezier 
polynomials allow us to easily impose both the monotonicity requirement and the 
matching of inflow and outflow cross-sectional areas. Consider 


B = {A € C\ 0, 1) |A(x) * AiMx) + qBiM + A wlt B 2 . J (x); 

x€ (0,1), gelAh,, Am)}, (35) 

where 2?,> is defined in (21). Thus, B is a one parameter set of curves in A. We 
restrict our optimization problem to this set B. 

Our final step in the approximation of Problem 2.3 is replacing the integral by a 
quadrature rule, with the set of quadrature weights and points {(c,,x, )}?„,. We now 
state the approximate design problem. 
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Problem 2.4 (Approximate Duct Design) Let {«,}?„, represent data for a de- 
sired transonic flow profile in the duet. We assume that the data and the approximate 
solution are given at the quadrature points, otherwise interpolation must be used. De- 
fine the objective function 

J, K lA) = t, «.[<(■*)-«,]’ (36) 

iml 

where u N (A) is an approximate solution to (26) with the cross-sectional area A. The 
approximate design problem is to find an A, € B such that 

J, N (A.) < J?{A) for all A € B. (37) 

Note that we can identify any A € B with the parameter q € Q s [Ak, A**] which 
uniquely represents it. Thus we can equivalently state the problem as to find q.6Q 
such that 

< J?(q) for all q € Q. (38) 

3 Sensitivity Equation Method 
3.1 Trust-Region Algorithms 

We shall use a trust-region algorithm for the optimization loop. The reason for 
selecting this type of scheme will be clear when we discuss the convergence properties 
in Section 4. This is a well known algorithm. However, we give a brief description 
below in order to prepare for the formulation of the sensitivity equation method. 

The quasi- Newton optimization algorithm produces a sequence of iterates which 
are obtained by minimizing a local quadratic model of the objective function. This 
model is constructed using the evaluation of the objective function J*{q k )-. its gra- 
dient V (q k ) and a secant approximation to its Hessian, H k at the current iterate 
The minimization of this model produces the next iterate q k +i, i.e. 

m k (q h+1 ) * min m*(q* + s k ) = min ( J g N (q k ) + VJ?(q k ) T s k + \s T k H k s^ . (39) 

Thus the next step is 

4*+i ~q k - H;'VJ?{q k ). 

It is well known that for sufficiently close initial guesses (and assumptions on the 
objective function), the iterates converge superlinearly to the minimum, q,. 
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However, the initial guess may not be in this superlinear region. Thus globaliza- 
tion strategies are employed to bring the iterates into the superlinear region. It is 
desirable to choose strategies which reduce to the quasi-Newton algorithm close to 
the minimum. One such strategy is a trust-region algorithm. In this algorit hm , a 
quantity 6, known as the trust-region radius, is used to measure the region in which 
the local quadratic mode), m*, is "trusted ’ 1 as an approximation of the actual ob- 
jective function, . Thus, the next iterate, q k +i> is now found by minimizing the 
model in this region, i.e. 


m *( 9 fc+i ) = min^ m k (q k + s k ). (40) 

where S k is the trust-region radius at the fcth iteration. 

A heuristic for changing the trust-region radius needs to be developed which in- 
creases S k when the model prediction is good and decreases S k when the model pre- 
diction is poor. One such strategy uses the ratio, 


(41) 


_ - 4TW) 

**kfafc) - m*(o+i) 

which is the ratio of the computed reduction to the reduction predicted by the model. 
If this ratio is small (or negative), then the model did a poor job of predicting 
and the trust-region is decreased. Whereas, if the ratio is near 1 , then the model did 
very well at predicting J** and the trust-region radius is increased. 

We present the resulting trust-region algorithm below. 


Algorithm S.l (Trust-Region) 

Select an initial guess <70 € Q, an initial trust-region radius So and constants 0 <>71 < 
T )2 < 1 and 0 < 7 ! < 1 < 7 ,. Compute J t N {q 0 ), VJ g N (q 0 ) and select or initialize H 0 . 
Do fc * 0,1,..., until “convergence” 


1 . Determine the approximate solution a k to equation (40). We chose the optimally 
constrained hook-step method [ 1 1 ] to do this. 

2 . If p k < T71, then set 6 ( 0 , 7 i£*) and q h+l m q k , jf(q k+} ) - J a N {q k ), 
VJfWk+x) * VJ*(q k ) a’id H k+1 - H„. 

3. If f)i <Pk< Vi, then set J fc+1 6 (0, 5*1 and q k+1 * q k + a*. Compute J*(q Hl ), 
VJ/W) and the update /f* +1 . 
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4 . If < pk, then set £*+i € [£*,72^*] and 9 *+i = 9 * + a*. Compute Jf*{qk+ 1), 
VJ]^(9* + i) and the update Hk+ 1 . 

Continue 


3.2 Design Sensitivities 


Is order to apply a gradient based optimization algorithm, such as the trust-region 
algorithm described above, we need to consider methods for computing the gradient 
of ■ In this discussion, we consider finding the gradient of (or a suitable 
approximation) with respect to the single design parameter q. This discussion can be 
easily extended to find the gradient of J* 3 with respect to multiple design parameters. 
A straight forward approach is to use a finite difference approximation, e.g. 



+ A<) — J?(,) 

A? 


(42) 


Unfortunately, this approach may not be practical for problems where the approx- 
imation of the PDE is computationally expensive, and is overly complex in shape 
optimization problems due to the necessity of computing mesh sensitivities. One way 
of alleviating the computational burden is to use design sensitivities, quantities which 
describe the influence of the design variables on the flow variables. For example, we 
can directly compute the gradient by differentiating (36) as 



(43) 


The quantity j^u N * is the design sensitivity for the discretized flow u N . 

There are several ways to compute this sensitivity. As above, one might use finite 
differences, yielding the approximation 

d S/ . g-f Aq) - u s (xnq) , „ 

ej u (*« «) * — ■ ' 2, — L±S1 - («> 

When the discretization is parameter dependent, it is easier to compute this approx- 
imation using, 


~u y (x,;q) « + £-M(gQAg;q 4 . A q) - tt^x/.q) 

Oq A q 

( 45 ) 
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in order to avoid interpolating back to the unperturbed mesh. This approach has the 
advantage that it may be possible to select a step size Aq using error estimates for 
«*. However, it is as computationally expensive as computing finite differences on 

A more efficient approach can be obtained by differentiating the simulation scheme 
used to approximate the flow (the discrete sensitivity approach). For example, in 
the FBS design problem, the simulation scheme (16) could be differentiated with 
respect to q , leading to a numerical scheme for terms like f^u N . Since the chain rule 
must be used to carry this out, the resulting scheme for the sensitivities contains 
terms similar to those found in the simulation scheme. Thus, the sensitivities ca n 
be computed efficiently along with the flow. A disadvantage of this approach is that 
when the discretization is parameter dependent, as in shape optimization probl em s, 
then derivatives of the discretization (terms like ^Af) need to be considered, see e.g. 
{ 20 ). 

An alternative approach is based on differentiating the original flow equation with 
respect to the design parameter and then approximating the resulting sen ntivity equa- 
tion- The result is (^v) , where the superscript JV refers to the approximation 

of the flow equation and the superscript M refers to the approximation of the sen- 
sitivity equation. Since this approach interchanges the order of differentiation and 
approximation, no mesh sensitivities are required. Furthermore, it has been shown 
(2] that applying the same approximation scheme to the sensitivity equation leads to 
similar computational advantages ss the discrete approach described above. More- 
over, additional computational savings could be obtained by applying a sch em e which 
takes advantage of the linearity of the sensitivity equation. A potential disadvantage 
of this approach, however, is that in general -fcu N qt ' M , even if the same 

approximation scheme is used for both the flow and sensitivity equations. 

However, if we consider the gradient of the infinite dimensional objective function, 

T q J{,) “ s /,’ M*i A > - 4 <*» Tg u{x ' A >- («) 

then using the sensitivity equation approach provides an approximation of this gra- 
dient, i.e. 

6 ( 8 \ SM • / v //» \ N ' U 

di J M * M “ 2 («*(*) - “’) (^ u ) («)• (47) 

Thus, we have reason to expect that this approach could produce feasible gradients 
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for the optimization scheme. These two sensitivity approaches are described in detail 
in later sections using concrete examples. 

3.3 Sensitivity Equation Method 

The sensitivity equation method couples a trust-region optimization algorithm with 
gradient evaluations provided by approximating the sensitivity equation. Thus we 
consider applying Algorithm 3.1 with the following quadratic model, 

M<lk+i) = H nuu^ Mik + s k ) = +$*** + H k s^j , (48) 

Note that we replace the quadratic model ms by to emphasize the fact that VJ* 
is approximated by p*, computed as {fcj} *"(?*)• 

The intent is to use the robustness of the trust-region optimization algorithm to 
compensate for the non-consistent gradients. The result is an optimal design method 
which is often more efficient and considerably easier to implement than current meth- 
ods. In the sections below, we discuss convergence issues and describe the implemen- 
tation of this method. 


4 Convergence Issues 

Definition 4.1 A numerical scheme is said to produce consistent derivatives with 
respect to approximations N (for the states ) and M (for the sensitivities ) if 

, NM 

(49) 


- (£/), <•> 


This is exactly the case for the discrete sensitivity approach, since one actually defines 
(computes) (^j)^ {■) to he 

Definition 4.2 A numerical scheme is said to produce asymptotically consistent 
derivatives with respect to approximations N (for the states) and M (for the sensi- 
tivities) if 




0, Vg € Go- 


(50) 


is satisfied as the approximations N and M art refined. 

We now consider the convergence of the sensitivity equation method. To begin 
with, we assume that the following hypotheses hold, 
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(HI) For a given qo in the design space Q, let Qo be an open convex subset containing 
the level set of J* at qo, i.e. 

A, = {» <= e | J/'U) < c Co S c. ( 51 ) 

(H 2 ) Jf* is bounded below 

(H3) J* 1 is Prechet differentiable on Qo 

(H4) The FVechet derivative of J* t denoted by VJ*, is Lipschitz continuous on Qo 
with Lipschitz constant L, i.e. 

l^/V) - V4V)|| 5*1 v,\ e e 0 . (52) 

(H5) The approximate gradient, y* is asymptotically consistent to VJ^ v (y*). 

(H 6 ) There exists a constant c\ € ( 0 , 1 ] such that 

< (-9k, »k) < IM!K|| Vfc as 1,2,. . . (53) 

(H7) There exist constants ca,c 3 € ( 0 , oo) such that 

-ct(d,d) < (H h d,d) <co(d t d) Vi * 1 , 2 ... ( 54 ) 


The following discussion parallels the proof given in (7} which treats the use of. 
trust-region algorithms with inexact gradient and function values. This discussion 
makes use of the fact that we seek the minimum of Jf* and have asymptotically 
consistent derivatives. 

Lemma 4.1 Under assumptions (H6) and (H7), Algorithm 3.1 produces iterates 
which satisfy 

Mlk) - tMfHi) > 5 C 1 M min jtf*, flSM j . ( 55 ) 
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(56) 


Proof Note that since ^*(?*) * J?{<lk), 

“ 0*(f*+l) = ~ ( 9k, »k ) ~ £ (HkSk,S k ) . 

Now, let it = ||«*||]^jj s a,d k , then o. solves 

a b*' <**> + \ a * (Hkdk,d h ) . (57) 

We can break this up into two cases, when (#*d*, dk) > 0 and when (//*<£*, d*) < 0. 
Case 1: Assume (JEM*, dk) > 0, then either 

fl _ <».*) 

’ ~ (Hkdk, dk) ’ 

in which case 


*<*)-*<*«) - <|^> M) ~ {Hidk ' dh) 

" 2{Hkdk t dk) -7T C3 
using hypotheses (H6) and (H7), or 


a.mS h 


in which case 


implies 


Sk< 


(g*»4) 

</***, d*> 


Mri-Mik+x) - -^(ff*,d*>-^|(H fc d*,dO 

- “^k {Skydk) + -6k {9k, dk) > ^Ci^UpfcH 

by hypothesis (H6). 

Case £: Assume d*) < 0, then o, = d*. Therefore 

0*(*k) - lM<?k + i) - ~«k(p*,d*)-^(^,d>) 

> Sk {9k, d k ) > c,6*)|p*U > ^c,^||p*||. 


A 
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(58) 


Lemma 4.2 Assume (HI) holds, then 

lim iof Hull > 0 and lim <5* = 0 

k-+ 00 


imply 

lim iHlUl - 1 

Proof It was shown {11] that, if || 4 *|| * 6k, then the solution 
s(pk), where 


(59) 

to (48) is given by 


•00 *«-(** + *!)“■* 

and in, is the unique zeal number that satisfies ||e(/x 4 )|| «= 6 k . Therefore, if 6 k 0, 

then nu -* oo (since Hi, is bounded, by (H7)). Thus ** -* “Ms’y*- A 


Lemma 4.3 Let J g N satisfy (HS), (Hf) and (H7), then the iterates satisfy 

I'M?*) - ^(?«i))-[j,"(?») - J/'tiMiJ < \ (c + L) ||.,|| > -( S1 - VJ„" (,*),.») . 

( 60 ) 

Proof Using the Cauchy-Schwartz inequality and (H3), we obtain 
(?*+]) -J*{qk) » (VJ*(g. + A»i),n)iiA 

- «*> + £ {^"(ft + A..) - VJ»(,.),..) iX 

£ n ) + £ | v^"(„ + u t ) - vj,"(„)| 

By the Lipschitz hypothesis (H4), 

< (v j, "(?.), ..) + £ iha», n y.„!i<fA 

- {vj,''(„),.,) + |i||n|| 1 . 

Thus, using (HT), 

iM<ik) - tMf*+i)] - [Jfiih) - J/'to+i)] 

£ - if,,;) - j .),«,) + jtll.,11 1 

£ “(y* “ '^# v (gi),..)+ ^ (c. + 1) ((.fcfl 1 
which completes the proof. A 
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Lemma 4.4 Assume satisfies (US), (HS) and (H4), and assume (HI) holds, 
then VJ** is bounded on Cq. 


Proof Let c be a constant such that J?(q) > c,Vq € Co (as guaranteed by (H2)). 
Assume to the contrary that there exists a point q € £o such that 

Define s - V where we choose a small enough so that $ + i € Qo- Then 

= - [ (VJ’/’W), *) <u - £ + as) - vj/'ff), j) i\ 

> - jiiur 

s f- L (1 ■ - § ) > £ (i - § ) [si wm - c)] 

This is positive for a € (0,2), thus J*(q) > J?{q + S), which implies f +*€£*. In 
addition 

J’/'W) - J’,"W+ ») > J, n M - C 

holds, but this is a contradiction dnce $ and $ + i are in Cq. A 


The <rem 4.1 Assume satisfies (HS), (HS) and (H4). Furthermore, assume the 
approximate gradient satisfies conditions (H5) and (H6) and that the update is con- 
structed so that (H7) holds. Then, for a sufficiently fine discretization, the sensitivity 
equation method produces a sequence of iterates such that 


UjnmfM *°. (61) 

Proof Assume to the contrary that liminf*^ ||y*j| > 0 and define $ k such that 

co s($ k ) * 

w MINI 

and v>k € Q such that 

f 0 sin(0*) = 0 

l a^(lfti + C08 ^)ifti) •*“<*») 1*0 ' 

Then (p*, u>*) ** 0 by construction, and 

(ft - vj/'UO.u.*) . . 

10 


If sin(0*) jt 0, then ||u>*|| * 1 and 

*k * IK II co 8 ( 0 *)jj~j| + sin( 0 *)u^ . (62) 

Let K denote the set of successful iterations, then 

Mih)-MQh+ 1) >m 
for each k € K. Lemma 4.1 implies 

3?(u) - J/W > I»in {«., 2M j . 

Since J* is bounded below, by (H2), the above condition implies lim*_ 0O ,* € jc £* — 0. 
Therefore, as £* is decreased in unsuccessful iterations, lim*-,^ 6* = 0. We now have 
the conditions for Lemma 4.2, and 

MINI ~ 

Thus cos(0*) » 1 and lim*-,^ sin(0*) a 0. 

Consider the expression 

1 _ = ~ ~ [Jfilk) ~ J?{<lk +\ )) 

Pk ” Mlk) ~ 'i’h(<lk+ 1 ) ’ 

by Lemma 4.3 and the definition of we get 

_ K°i + jjlkf - U - vj?(ik),>k) 
n ~ (-«».*») - J ' 


Using hypothesis (H7), 


(-0*1 9 k ) 

Substituting expression (62) and using ||s k || < 0*, we get 

1 ~ ( 0 * ~ v >7/ r (0»)»0*) ~ ( 0 * - VJ a N (<Ih),Wk) sin(0») 

«0*l|cos0* 

By Lemma 4.4 and the Cauchy-Schwarz inequality, (yj t N (q k ), w*} is bounded and 
we consider the limit as k -* 00, 

1 . , t , to i * - $>>■*) £ 1** - 

*-«• *-<» ||0fc|| 3 IK || 
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Since lira inf *_<» ||pi,|| > 0 and <7* is asymptotically consistent, we can select a suffi- 
ciently fine discretization such that 

lkl-/i*<l- 173. 

k—oo 

Hence, pk> which implies d*+i > a contradiction. A 


5 Duct Design Problem 

In this section, we use the duct design problem to illustrate the implementation of the 
sensitivity equation method. To begin with, we will introduce the discrete approach 
for finding design sensitivities in order to compare it with the sensitivity equation 
approach. 

5.1 Discrete Sensitivities 

To obtain an algorithm for the sensitivities ^u' v ($) as ($)}* , the system of 
nonlinear equations ( 32 ) is differentiated, yielding 

=0. < 63 , 

where P ,+j /j is determined by the scheme used to compute the flow. If the Enq uis t- 
Osher scheme was used, 



f « 

/i+i 


F% n - • 

h 

0 uj' < U, < tijii* 

. h + fi+l «J>1 < U. < of, 

( 64 ) 

or if the artificial viscosity scheme was used, 


*VX /3 * \ (/>+> + /i - * ( 1 - 

l) . («) 

where /,• = /(«?,£«?), 



'(■ 


(66) 

and 



1 ( u % u ' A 'i A ) - i (uc 4 ) (’■-*) + (\h A ) 

( ,+ £)r- < 67 > 
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This differentiated scheme can now be used to compute 
5.2 Sensitivity Equation 

We now present the implementation for the sensitivity equation approach. We begin 
by differentiating the flow equation (26) with respect to the parameter q. Thus 


s^(“4“) +# ( o, I“’ a 4 x ) =0 


b 


and 


r<»-° 




(69) 


is the sensitivity equation for this problem. Note that the sensitivity equation is a 
linear equation with variable coefficients (determined by u). There has been little 
analysis of numerical schemes to approximate equations of this type. However, for 
this two point boundary value problem, the same numerical schemes (Enquist-Osher 
and artificial viscosity) provide convergent algorithms. As in the approximation of 
(26), we consider U“); to be the average sensitivity in the ;th cell. A system of 

nonlinear equations for (^u) "(«) »= {(£“)" w}. can be found by integrating 
//M *' over each cell, * 

/ (u(*, 4 $), £u(gj 4 j)) - / (ufo - &), &t ;( x, - &)) 


+ 1 


” 0 ’ ( 70 ) 


j = 1, . . . , A r , where we assume A and f^A are nearly constant over each cell. As 
before, the terms / (u(xj 4 $), ^t*(«> + £)) are replaced by the cell center values Jj 
and Jj+\. Using the Enquist*Osher scheme, we obtain 


pBO 

r i + i/s 


f /*. 
h 

/(«#> (£«)#) 

fi + />+ 1 + /(«„ (£«),) uj+ t <u,< uf, 


> V, 


uf < u t < u 


,s . 

‘>+1* 


(H) 


and obtain 


fiAV 


H /w,+/i ' a (S ,, L'K u )!}) 


(72) 





for the artificial viicoeity scheme. It is obvious that the approximation of the sensi- 
tivity equations depends on the approximation of the flow equations. As described 
earlier, we use the notation to represent using scheme N to approximate 

the flow equation and scheme M to approximate the sensitivity equation. 

5.S Convergence Results 

The convergence result provided ir Theorem 4.1 can be proved for the case when 
the artificial viscosity scheme is used to approximate the flow and the Enquist-Osher 
scheme is used to approximate the sensitivities in Algorithm 3.1. For this problem, 
we assume (the (Hi) in Theorem 4.1) that 

C “ [Ain, Aeut] , Qo = (Ain; Amu) , 
and 

£. = {«« a | j?" («) < } . 

The objective function J^ AV given above is obviously bounded below (by zero if all 
of the quadrature weights are nonnegative) satisfying (H2). The hypothesis (H3), the 
differentiability of 

- £<* («*"<*,«) - so *)) 1 

1*1 

on Co and hypothesis (H4), the Lipschitz continuity of the derivative, follow from the 
following 

Lemma S.l The approximate solution u Nav is differentiable and the derivative is 
Lipschitz continuous on Qo. 

Proof The approximate solution, u N * v is the root of the nonlinear equations 

W , q) * [F# /3 (u^v , ,) - F% /7 (u^v , ,)] + g . ( u ^v , ^ o, (73) 

where Fff x/ a and g are C* functions of their arguments (for u Nav > 0). Then by the 
implicit function theorem, the map 

q-*vP* v (q) 

is Lipschitz continuously differentiable. A 

We point out that the differentiability of the approximate objective functional is 
strongly dependent on the discretization scheme used in the approximation. For 
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example, the objective functional associated with a Godunov approximation of the 
flow is not differentiable, a result of matching a parameter dependent discontinuity on 
a discrete set of points [4]. Finding feasible optimisation strategies for this problem 
has been the focus of recent work, see e.g. |4], (19] and (23). However, for the 
purpose of this discussion, the artificial viscosity scheme provides a smooth enough 
approximate objective function. 

The hypothesis (H$) is guaranteed (for some discretization level) by the asymptotic 
consistency shown below. 

Theorem 5.1 For the one dimensional Euler equations, the derivative 
where the flow is approximated using the artificial viscosity approximation and the 
sensitivities are approximated using the Enquist-Osher scheme, is asymptotically con- 
sistent to -§^J * AV . 

Proof Consider the norm used in the definition of asymptotic consistency: 



The first term on the nght hand side vanishes since using the artificial viscosity 
scheme for approximating both the flow and sensitivity equations leads to consistent 
derivatives. The last two terms go to zero as the approximations JV^v, Mav and M^o 
ate refined, since the artificial viscosity and Enquist-Osher schemes converge when 
used to approximate the sensitivity equation, is the exact solution to 

the sensitivity equation given u Nav . A 

The hypothesis (H6) can be enforced by the optimization algori thm by rejecting 
steps which violate this condition and shrinking the trust-region radius. This proce- 
dure eventually creates a step which satisfies (H6), since the limit of this procedure 
would produce a step in the steepest descent direction. 

Finally, (H7) can be enforced by the secant update strategy. Therefore, we have 
shown that these approximation schemes satisfy the conditions of Theorem 4.1. Nu- 
merical computations using these sensitivity schemes are provided below. 
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Figure 3: Deeign Seoaitivity Approximations Using Enquist-Osher Scheme 

5.4 Numerical Results 

The sensitivity of the velocity with respect to the Bezier parameter, q, is presented us- 
ing the numerical schemes described above. For this computation, the cross-sectional 
area corresponds to an element of B (see (35)) with q ~ 1.37125. The interval [0, 1] 
is divided into 45 cells. In Figure 3, the sensitivity solution using the Enquist-Osher 
scheme to compute both the flow u N *° and the sensitivity (fo)* 0 '”’ 0 b compared 
with the closed form sensitivity solution. In addition, the sensitivities computed 
via finite differences of Enquist-Osher solutions using a finite difference step size of 
A? * (1 x 10"*) q are also provided. Excellent agreement is seen for both of these 
methods. The only discrepancy is in the cell to the left of the shock, where numerical 
dissipation appears in the flow solution. 

The corresponding design sensitivities which are computed using only the artificial 
viscosity schemes are shown in Figure 4. As above, the agreement is excellent except 
where dissipation errors appear in the flow approximations. In this case, these errors 
appear over more cells near the shock. 

Note that the computation of these sensitivities were performed efficiently, rela- 


25 





Figure 4: Design Sensitivity Approximations Using Artificial Viscosity Scheme 


tive to the cost of & flow approximation. The flow approximation requires solving 
a system of nonlinear equations. The sensitivity approximation, on the other hand, 
only requires solving a linear system since the sensitivity appears only linearly in the 
definition of / and g. Moreover, if the Newton method is used to solve the nonlin- 
ear system, then the linear system is already available in factored form. Therefore, 
the sensitivities can be computed using less computational time than required for 
one Newton step. Computational efficiencies such as this can be missed if the flow 
algorithm is simply differentiated. 

Note that as long as is bounded, 


/ ( u *’U“l) “('-I) 


since fi » t»J. Thus, one observes that the numerical algorithms to compute ei- 
ther t-y*o ox (Jju) *°' *° are equivalent. This leads to the fact that using the 
Enquist-Osher scheme to approximate both the flow and sensitivity equations pro- 
duces consistent gradients. In addition, it it easily seen that using the artificial 
viscosity scheme to approximate both equations also produces consistent gradients. 



Table I; A Comparison of Gradients at the Optimum for Various Mesh Sizes 


N 

(• 

mmm 


15 


I'MIjrftl 

WMiEmmat 

45 

1.3437 


•0.001566 

135 

tml 

0.002485 

•0.000012 

225 

ttsa 

0.002478 

0.007802 

315 

1.3553 

0.002645 

0.026731 

405 

1.3554 

0.002816 

0.001584 


However, if the artificial viscosity scheme is used to approximate the flow and the 
Enquist- Osber scheme is used to approximate the sensitivity equations, the gradients 
are not consistent but asymptotically consistent. 

Numerical results for this asymptotically consistent case are provided in Table I. 


6 Forebcdy Simulator Design Problem 

We now describe the implementation of the sensitivity equation method for the fore- 
body simulator design problem described in Section 2. As in the duct design problem, 
we begin by presenting the equations which comprise the discrete sensitivity scheme 
in order to compare and contrast the two methods. Unlike the duct problem, we 
have no theoretical convergence results for the FBS design problem. However, the 
numerical experiments below show that the SEM still converges. 


6.1 Discrete Sensitivities 


Differentiating the numerical scheme (16) with respect to a design parameter, repre- 
sented by q , leads to the following scheme: 


I + AM, .4" - V,^ 1 -r »{•>)£,/„] x 


[J + - ?,(*?> + a!<5” = 

[/ + AM.fi* - V,(*f + *<*>)A,^] A$* 


- [/ + AM, A* - + *<*>)A, J M \ x 

f AM,lfi* - V,(±*m + 1*W)A V J M - V,(«f» + 


a& 
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+ AtV<(JU< ! > - 1 *<*>A ( V < )A < (J a ,(5”) 

+ A«V,(*f - 

+ A<V,(l*f - ^*<*>A,V,)A,(^(5") 

+ A(V,(*B) - *<«A,V„)A,1(J m <5*). (74) 

The equation representing the boundary conditions axe also differentiated. Note 
that the above sensitivity scheme requires derivatives of the mapping, (denoted 
as mesh sensitivities) and the dissipation terms, and ^9 (4) . Evaluation of 

TjM is given by differentiating the scheme which determines M , see e.g. [20]. Other 
methods for approximating j- q M have also been investigated, see e.g. [25]. We see 
from (74) that terms containing these expressions represent a signifi cant portion of 
the computational effort, aside from the fact that and themselves 

need to be determined. 


6.2 Sensitivity Equation 


The sensitivity equation approach to computing design sensitivities is presented be- 
low. To begin with, we differentiate the Euler equations and associated boundary 
conditions with respect to the design parameter q, which leads to: 


where 


df ; 8G t 
9x + 9y 


— 0 




£ 

Bq 


uQ + uQ t + 


0 


l h Pu + p if‘ \ 


d 

G ' ■ Tq'Q + + 


• 1 
0 


' *• 1 

0 

O sx 



, Vf a 

£(?») 



** J 


(75) 
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and where 



9. . 6 pu 

dq( pU ) “ dq P p J 


!(> 




since pjLQ. 

We are now free to apply any appropriate scheme to solve (75). In particular, it is 
possible to use a method which takes advantage of the linearity of the sensitivity equa- 
tion. However, in this work, the same scheme used to solve the flow equations is used 
to approximate the sensitivity equations, which leads to an efficient computational 
scheme as in the discrete approach (2). This scheme is described below. 

This equation may now be transformed to generalized coordinates, so that the finite 
differencing can be done more easily. It makes sense to use the same transformation 
(which is equivalent to using the same mesh) that was used in the solution of the 
Euler equations. Thus the resulting system is 




o, 


where 



» « 
0 


’ 0 ' 

f. « UQ, + U,Q + ~PJ^i 

s 

g 

+ PJm‘ 

0 

0 


u 


.V 


■ 0 1 


9 « 

0 


fe 

g 

+ PJ2 

0 

0 


V 

! 

. v < . 


where 

U = V( (u,v) and U t ■ Vf • 
and 

V * Vij • (ti,v) and V t « Vfj* 
It can be shown that 


(76) 


A 


dF dP. . b 

55 = 5^, »<t b 


dO ao, 

dQ~W,' 
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I 

» 

I 

I 


so th&t ths discretization has the same factored form as the Euler equations, thus 
[/ + - V,(*f> + *J 4 >)A t y*<] X 

[/ + Atf,S" - V„(*<’> + *<<>)A,J m ] A(l<})* = 

- m t (F,r - Ats,(6,r 
+ AlV,(*fl - *<«A < V < )A 4J M Qr 
+ - *< 4 >A,V,)4,M«0)* (77) 

Since the left hand side matrices are the same, a right hand side vector needs to 
be formed for each design sensitivity. In addition, the boundary condition type is 
the same for both the Euler and sensitivity equations. The boundary conditions are 
determined using implicit differentiation. 

Note that this scheme is similar to the discrete sensitivity approach. However, since 
the approximation is applied after the differentiation, there are no mesh sensitivity 
or dissipation sensitivity terms. The other obvious difference is that the boundary 
condition on the parameter dependent boundary is different. 

6.8 Boundary Conditions 

The boundary conditions for the sensitivity equation (75) are provided below for the 
case where the forebody simulator is described by a two parameter Bezier curve (18)- 
(20). Extensions to other forebody descriptions will be obvious. The appropriate 
conditions are obtained by differentiating the corresponding boundary conditions for 
the Euler equations. For example, at the inlet, we flow Q* is prescribed and will not 
vary as the forebody parameters q * (g l ,g a ) axe changed, thus 

Qi “0 

at the test cell inflow. The walls are treated in a similar fashion. However, the 
boundary condition at the forebody simulator surface requires more attention. This 
is because the points where the condition is evaluated are parameter dependent. 

We study the treatment of condition (5) in detail. The normal vector to the 
forebody surface is 
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Thus, the boundary condition (5) can be written as 

— n(r*(s),r„(s;g);g)^r if (e;g)"f v(r*(s),r|,(e;g);g)^r*(a) * 0. (79) 

The corresponding sensitivity equation boundary condition for the first parameter, 
g 1 , can be obtained via differentiation, i.e., 

~Jq* v ( r *( a )> r *( e; f ) + ( r *(«; «)> r *(«); 9) * 

(F*( g )i r*(s» q)\ ?} J»(s; ?)g^r *(*; q) 

+ «( r.W,r,(s;g);0^ i r,(.;f) 

“ (M*). r,(a; g); g) ~r v (s; fl)^r,(s). 

This is simply a nonhomogeneous version of condition (5), namely, 


( d d \ . d d . d 




Using the same techniques, the boundary conditions corresponding to (6) are: 

a 5 a_ $_ A a a 3 _ a 3 a _ a_ 
dxdy Ut dq l v d» v + a* u ‘asag> r> ~ a? tt ‘a? r ^ r# 

a / a \ a 3 a_a_^a a 3 _ a 3 a^a,, 

an lag 1 v ~ a^av^dg^as^ * dx Vt dsdq' r * 


The analogous boundary conditions for g 3 are obvious. 


(80) 


6.4 Numerical Results 

The sensitivity equation approach, which computes design sensitivities for the two 
dimensional Euler equation is illustrated below. In this implementation, a right hand 
side vector for each design sensitivity is formed along with the corresponding vector 
for the flow approxi m ations . The updates for the flow and sensitivity variables are 
obtained simultaneously, exploiting the fact that the left hand side matrices are the 
same. 

The design sensitivities with respect to the first Bezier parameter g 1 were computed 
for a forebody described by the curve 

f-(t(j),«j)). *6(0,11, 
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where 


x(a) - 0.05e I »(e) + 0.1«B,3(<) + 0.55fl, l3 (i)+1.0F9^(e), 

v( 8 ) ** r 9 Bo*(s) + + g 2 £ 3 i 3 (a) + rfcBj4(j), 

g 1 = 0.1, f 2 = 0.15, T, = 0 and Fi = 0.2. This curve is twice as long in the z-direction 
as the admissible forebody simulators given in B (see (18)). Under a uniform inlet 
flow profile described by the inlet Mach number, M, = 2.0, the approximate flow 
variables and sensitivities are computed on a 43 x 49 mesh. The sensitivity of the 
e-component of momentum with respect to the Bezier parameter computed using 
the sensitivity equation approach and the finite difference approach (for 4 different 
step sizes) are plotted along the outflow plane in Figure 5. The corresponding plots 
for the Energy sensitivity are provided in Figure 6. Observe that the step size of 
0.00001 produces noisy sensitivity values close to the forebody (presumably due to 
round-off errors). A larger step size of 0.01 gives the best results (when compared 
to the sensitivity equation approach) near the shock location. The best qualitative 
behavior appears when the step size is 0.001. These figures demonstrate the difficulty 
of obtaining a satisfactory step size at all resolution levels in the flow domain. 

A model forebody simulator design problem is discussed below. To begin witn, 
we seek the optimum value of the inlet Mach number and two Bezier parameters 
( (9 l >« 3 ). describing a shortened forebody simulator in the admissible set B) which 
minimize the approximate cost functional J” (given in equation (25)). The flow 
data 0 to be matched is given by the flow Q N corresponding to the forebody shape 
f described above. We point out that the artificial dissipation in the flow solver 
produces a “smearing” effect on the flow variables. Therefore, based on the results 
for the duct design problem, we expect a sufficiently smooth approximate coat func- 
tional. Furthermore, the comparison of the sensitivities in Figures 5 and 6 lead us to 
believe that the sensitivity equation approach may produce asymptotically consistent 
derivatives. 

The sensitivity equation method was applied to the FBS design problem with 
initial values of the parameters: M a = 2.0, q l * 0. 10 and q 2 * 0. 15. These parameters 
correspond to those used to generate 0 (even though that forebody is longer). We 
present the iteration history in Thble II. Observe that there is a drastic reduction 
in the approximate cost functional in the first three iterations. The iteration history 
for the i- component of momentum is given in Figure 7. Note that the front end of 
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Figure 5: Compariaon of Sensitivitie* at Outflow: x* Component of Momentum 
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Figun 6: Comparison of Sraoitivitie* at Outflow: Energy 


93 



Table 11: Shortened Forebody Optimisation 


| Iteration 


■ran 

■nm 

msmirmarm 

iszssa ] 


|M|J 




KciEUi 



0.14008 

rl 


11.6285 





04882 

3.7955 



0.80765 

Iff 7?Ti 

0.2334 



2.01027 

0.30189 

0.14007 

0.2806 

0.5963 



0.29367 

0.14787 

0.2289 

0.6861 


o 

0.28891 

0.15564 

0.2271 




0.29011 

0.15921 

0.2249 

0.1518 


Ji? 1 

0.29278 

0.15821 

0.2237 

^jTiTT7TM 



0.29420 

0.15689 

0.2232 


10 

2.01952 

ITOI 

0.15604 

0.2280 

0.0275 

11 

2.01994 

liiilll 

0.15608 

1 IvvlM 

0.0173 

12 

2.02006 

0.29415 

0.15609 

0.2229 

0.0158 


the fore body simulator becomes more blunt during the first two iterations in which 
a stagnation region is set up in front of the FBS. This has the effect of moving the 
shock forward, which comes close to the shock location created by the long forebody. 
The remaining iterations are used to “fine tune” the solution near the FBS. The 
comparison of the optimal forebody simulator to the flow generated by the long 
forebody is displayed in Figum 8. Notice that the shock location is the same in both 
flows. 

In the optimization above, the initial Hessian was computed using forward dif- 
ferences. This adds some initial expense in the hope for fewer iterations. However, 
without this technique, using the identity matrix as the initial Hessian, the iteration 
converged m fifteen iterations. Therefore, neither technique showed an advantage. 

6.6 Conclusions 

While no rigorous proof of asymptotically consistent gradients has been shown for 
Euler equations, numerical evidence in (3] suggests that the gradients may indeed 
be asymptotically consistent. Similar numerical evidence exists for finite element 
approximations of the Navier-Stokes equations [5]. 
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Final Converged Solution Original Long Forebody (data to be matched) 

Fifun 8: Companion of Optimal Solution with Data 
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